Finding the center reliably: robust patterns of developmental gene expression 
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We investigate a mechanism for the robust identification of the center of a developing biological 
system. We assume the existence of two morphogen gradients, an activator emanating from the 
anterior, and a co-repressor from the posterior. The co-repressor inhibits the action of the activator 
in switching on target genes. We apply this system to Drosophila embryos, where we predict the 
existence of a hitherto undetected posterior co-repressor. Using mathematical modelling, we show 
that a symmetric activator-co-repressor model can quantitatively explain the precise mid-embryo 
expression boundary of the hunchback gene, and the scaling of this pattern with embryo size. 
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During embryonic development, cells have to differenti- 
ate in a manner that is dictated by their positions within 
the developing embryo. Turing conjectured that the posi- 
tional information for differentiating cells is provided by 
molecules called morphogens that self-organize into spa- 
tial patterns 0|. Today a number of proteins have been 
identified that act as morphogens. However, in contrast 
to Turing's original conjecture, morphogens are usually 
produced at a localized source, after which they diffuse 
into the surrounding tissue. Here, morphogens are also 
degraded and a concentration gradient is thus formed. 
The differentiating cells respond via patterns of gene ex- 
pression that depend on the morphogen concentration . 
However, the biochemical processes that give rise to the 
morphogen gradient are prone to variations: the rates 
of protein synthesis/degradation, for example, will vary 
from embryo to embryo. Variation in these parameters 
might modify the morphogen gradient and thereby in- 
duce error in the positional information transmitted to 
the differentiating cells. Nevertheless, cell differentiation 
is often exceedingly precise Q. How embryonic devel- 
opment is robust against fluctuations in the underlying 
biochemical processes is still poorly understood. 

Recently, several mechanisms have been proposed to 
explain robust embryonic development 0, IE @- They 
rely on the insensitivity of the morphogen gradient to 
embryo-to-embryo variations in the morphogen synthe- 
sis rates 0, IE El- However, these mechanisms cannot 
explain recent experimental observations on the Bicoid- 
Hunchback system in the developing Drosophila em- 
bryo bicoid (bed) messenger RNA (mRNA) is de- 
posited at the anterior pole of the embryo, where it 
forms a localized source of Bed protein. The Bed pro- 
tein diffuses away from the pole and is degraded, thereby 
forming a concentration gradient (see Fig. la). A sim- 
ple 'threshold' model would postulate that downstream 
genes such as hunchback (hb) are activated only when 
the concentration of Bed is above a certain threshold Q • 



However, while the hb expression domain is highly pre- 
cise, the Bed morphogen gradient exhibits large embryo- 
to-embryo variations [E Q- This rules out the above 
mechanisms, which rely on the robust formation of the 
morphogen gradient itself. Moreover, experiments reveal 
that the Hb boundary, rrnb > scales with the embryo length 
(EL): x Hh = 0.49 (±0.01) EL 0. Such scaling of expres- 
sion domains cannot be explained by a single morphogen- 
gradient model, in which positional information is only 
emanating from one end. Indeed, this observation sug- 
gests that the formation of the hb expression domain is 
regulated by two or more morphogens, sourced at oppo- 
site ends of the embryo. Earlier work of Wolpert Q and 
Mcinhardt 8] already suggested that such models might 
provide size regulation. 

Motivated by recent experimental observations 

HE3, 

which show that the activity of Bed itself can be reduced 
by interactions with other proteins called co-repressors, 
we predict the existence of a second, as yet unidenti- 
fied, morphogen that diffuses from a posterior-localized 
mRNA source. The co-repressor binds to Bed, thereby 
inhibiting its ability to activate hb expression in the pos- 
terior half of the embryo (see Fig.QJi). We show that such 
a symmetric activator-co-repressor module can locate the 
embryo center reliably. Our model thus generates ro- 
bust patterns of gene expression that scale with size, a 
result of potentially wide importance. Our calculations 
reveal that the module is highly robust against correlated 
embryo-to-embryo variations in the synthesis rates of the 
anterior activator (Bed) and posterior co-repressor. Fur- 
thermore, variations in the synthesis rate of Bed that are 
uncorrelated with that of the co-repressor, are also fil- 
tered far more than in a single-morphogen model. Iden- 
tifying the proposed posterior morphogen may require a 
large-scale search for maternal mRNA localized to the 
posterior pole, which does not correspond to an already 
characterized maternal gene. 

Mathematical Model. The outline of our model is 
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shown in Fig. ^ Bed and co-repressor protein are in- 
jected at rates J B , J b at opposite ends of the embryo, 
where at each time instant the injection is concentrated 
at a randomly chosen location x B (t), x B (t) between 5/im 
and 30/im from the respective poles. The newly trans- 
lated proteins diffuse away from the poles with diffusion 
constant D, and are degraded at a rate p.. Bed and the 
co-repressor bind with rate v to form a complex also with 
diffusion constant D. We assume that the reverse reac- 
tion occurs at a far slower rate, which we neglect. Only 
active Bed (unbound to co-repressor) binds to the hb con- 
trol region on the DNA to activate its transcription. Ex- 
periments suggest that this is a cooperative process |l ll , 
which we model with a Hill function with a coefficient 
of 3 |ll|. However, the precise form of the activation 
function is not very important: taking other functional 
forms or changing the Hill coefficient to unity, for in- 
stance, does not alter our results, as the concentration of 
active Bed is reduced to very low values at mid-embryo 
due to the co-repressor (see also below). Hb also auto- 
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FIG. 1: a: Sketch of Drosophila embryo showing localized 
Bcd/co-repressor sources at opposite poles, and co-repressor- 
Bcd binding, b-g: Positional information for Bed and Hb 
density profiles; data from 100 simulated embryos. Maxi- 
mum densities normalized to unity. Mean density profile and 
standard deviation for: b: Total Bed (Bed plus complex); c: 
Unbound Bed and co-repressor; d: Hb; in b)-d), 10 typical 
individual density profiles are also shown, e: Distribution of 
length scale A describing exponential decay of total Bed den- 
sity profile, f: Distribution of positions XBcd and xnb where 
each Bed and Hb density profile crosses given threshold, 0.17 
for total Bed, 0.5 for Hb. g: Positions XBcd and xuh as func- 
tion of embryo length (EL); Bed: circles, Hb: squares. Green 
straight line is xub = 0.49 x EL. 



activates its own synthesis via 'synergistic' interactions 
with (active) Bed 0, ^3 • mRNA has diffusion con- 
stant Dhb, and is spontaneously degraded with rate /Jhb- 
Hb protein is manufactured at a rate proportional to the 
mRNA concentration, has diffusion constant -Drib and is 
spontaneously degraded with rate /iHb- These processes 
give the following equations, where [B] is the unbound 
Bed density, [R] the co-repressor density, [C] the Bcd- 
co-repressor complex density, [hb] the hb mRNA density 
and [Hb] the density of Hb protein: 
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Here, we have used a simplified Id model, where x mea- 
sures the distance from the embryo's anterior pole. Zero 
flux boundary conditions were imposed at both embryo 
ends. The wild- type parameters used were (see also jjj): 
EL = 460/xni, D = 10^m 2 s _1 , D hb = 0.5//m 2 s"\ -D H b = 
1/imV 1 , J B = Jr. = 0.08s" 1 , n = 0.0007s" 1 , fi hb = 
0.004s" 1 , ^Hb = 0.006s" 1 , v = 3/mis" 1 , a = 0.006s" 1 , 
13 = 0.004 / um" 1 s" 1 , r) = 1.0, 7 = 0.01, K = 0.045 / um" 1 . 
We use the same diffusion constant/degradation rate for 
the complex as for unbound Bed; this ensures that the 
total Bed profile decays exponentially (as seen in experi- 
ment [3]), but is not important for the basic center find- 
ing properties of the model. Spatiotemporal fluctuations 
of the components within a single embryo, which could 
affect robustness were also modelled but found to 
be less important than embryo-to-embryo variability and 
were therefore neglected. To simulate the embryo-to- 
embryo variability, the embryo length fluctuated by 10% 
0. Other parameters were varied by imposing a Gaus- 
sian distribution about the above average values, using 
the following standard deviation divided by the mean: 
J B : 0.05, J R : 0.05, fi : 0.3, fx hb : 0.3, p Hh : 0.3, a : 0.3, : 
0.3, K : 0.3. The variations in the protein degradation 
rates were assumed to be correlated, while we studied 
both correlated and uncorrelated fluctuations in the Bed 
and co-repressor source strengths. Initially, the system 
was empty of proteins and hb mRNA. The system reaches 
a steady state after ~ 1 hour, well before the window in 
which the relevant experiments were performed yMTj. 

It is conceivable that an additional activator dif- 
fusing from the anterior competes with the co-repressor 
for Bed binding. This, however, does not significantly 
change the results. Models in which the co-repressor 
and/or the Bcd-co-repressor complex can also bind to 
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FIG. 2: a: hb mRNA profiles in simulated Bcd(A52-56) 
mutant (y — 0.003/ims" 1 ). Inset: mean/standard devi- 
ation for position of hb boundary as function of Bcd/co- 
repressor binding rate v (closed symbol denotes wild type), b: 
Mean/standard deviation for Hb profile in mutant with mal- 
functioning Hb protein, c: Mean/standard deviation for Hb 
density profile in simulated staufen mutant (see text). Data 
from 100 simulated embryos; maximum densities normalized 
to unity; 10 typical individual profiles are also shown. 



the DNA, but cannot then initiate hb transcription, again 
produce very similar results. Schemes in which Bed and 
(co-) repressor only interact on the DNA, for instance by 
mutually excluding each other's DNA binding, provide 
scaling, but are less robust to variations in their param- 
eters than the model proposed here. Furthermore, we 
have neglected the posterior morphogen nanos because 
if nanos is removed, the key properties of accuracy and 
scaling remain |3j. Lastly, while the co-repressor inacti- 
vates Bed at mid-embryo with respect to hb expression, it 
does not impair Bed's capacity to activate other genes in 
the embryo's posterior, such as giant and knirps [15J. In- 
deed, the Bed domain that interacts with the co-repressor 
plays an inhibitory role in hb expression, but a positive 
role in knirps expression |l6| . 

An alternative model has been suggested 0]; ^ as- 
sumes that the mRNA-localization protein Staufen trans- 



ports hb mRNA away from the poles via a counter- 
intuitive mechanism relying on an increased diffusion rate 
of hb mRNA upon binding to Staufen. However, scaling 
is found over only 1-2 minutes rather than the required 
hour or more. Moreover, it predicts a gradient of Staufen, 
which is in contradiction with experiment ^^|. Hence, as 
proposed here, a different mechanism may be needed. 

Results. Simulation results of the above model, with 
uncorrelated fluctuations in the source stengths, are 
shown in Fig.^ The density of Bed plus Bcd-co-repressor 
complex decays exponentially away from the anterior 
pole (Fig. \£}p) , in agreement with experiment Q . The 
total Bed profile shows considerable embryo-to-embryo 
variation, with the characteristic length scale A = y/D/fi 
of the profile (Fig. and the position where the density 
crosses a threshold (Fig. [If ) , both showing large fluctu- 
ations, similar to experiment 0,0- However, the pro- 
files of the co-repressor/unbound Bed fluctuate far less 
from embryo to embryo (Fig^); the fact that they bind, 
forces their unbound densities to very low values at mid- 
embryo. The active Bed is then able to precisely acti- 
vate hb: the hb mRNA (data not shown) and Hb pro- 
tein (Fig. njl) densities both pass their half-maximum 
values at xm> with little embryo-to-embryo variability, 
XHb = 0.49 (±0.02) EL, see also Fig. [Tf. The concentra- 
tion profiles switch between their highest and lowest val- 
ues over a length scale of about 0.1 EL. Both the average 
position and variation of the Hb boundary, as well as its 
steepness, are in agreement with experiment 0,0- I m_ 
portantly, the Hb boundary location scales with the em- 
bryo size, unlike the Bed profile 0,0 (Fig.Qfe). Further- 
more, scaling is achieved even before the protein/mRNA 
concentrations reach steady-state, adding additional ro- 
bustness to the model. Finally, scaling can be obtained 
even when the source strengths and A are not the same 
for Bed and co-repressor (data not shown). 

A critical test of our model comes from varying the 
bed gene dosage. When the dosage is halved (i.e. Bed 
injection rate halved), the Hb boundary moves to xub = 
0.39 (±0.02) EL. Simulations with twice and three times 
normal dosage, but assuming only 50% efficiency for the 
extra bed genes as suggested by experiments 0, (i-e. 
1.5 and 2 times the wild type Bed injection rate), move 
the boundary to aiHb = 0.56 (±0.025) EL and xnb — 
0.60 (±0.03) EL, respectively. Both the results on the 
mean and the variance agree with experiment In- 
terestingly, the accuracy of hb expression is only weakly 
affected by changes in the bed gene dosage: the standard 
deviation in xm> only increases by 50% when the Bed in- 
jection rate changes from 0.5 to 2 times the wild type Bed 
injection rate. More importantly, the average hb bound- 
ary moves by less than predicted by a simple threshold 
model : the co-repressor gradient makes the boundary 
positioning more robust to variations in the protein pro- 
duction rates. Thus, the posterior morphogen not only 
provides scaling, but it also plays an important role as a 
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buffer against gene expression fluctuations. 

Our model is also supported by several mutants. The 
Bcd(A52-56) mutant inhibits the ability of Bed to inter- 
act with co- repressors 0. We have modeled this by re- 
ducing the Bcd/co-repressor binding by a factor of 1000. 
Fig. 0i shows that this leads to a pronounced posterior 
shift in the hb expression pattern, in agreement with ex- 
periment We also predict that its precision will be 
lost: the hb boundary fluctuates strongly as the poste- 
rior co-repressor is no longer able to buffer against Bed 
fluctuations. Experiments have also been performed in 
mutants that make detectable but malfunctioning Hb . 
Normally Hb auto-activates its own synthesis [11 [3] i an d 
our model includes this effect. If, however, we assume 
that in this mutant Hb cannot stimulate its own synthe- 
sis (corresponding to rj = 0), then the Hb boundary shifts 
anteriorly to x H b = 0.46 (±0.02) EL (Fig.^j), while re- 
taining its scaling properties. These results all accord 
with experiment [3J, LLSj. Further support comes from 
staufen mutants. Staufen play s a central role in polar lo- 
calizing maternal mRNAs |l8l llfll ] . In the stau HL / stau r9 
mutants, the precision of the Hb boundary is lost and 
many of the Hb profiles are shifted anteriorly 0. In 
contrast, the Bed profile is not altered in the stau HL 
mutants Q . Taken together, these observations strongly 
suggest that in the stau H L / stau r9 mutants, co-repressor 
mRNA is no longer properly localized at the posterior 
pole. We have studied this scenario by injecting a ran- 
dom fraction of the co-repressor close to the posterior 
(as described above), with the remainder randomly dis- 
tributed throughout the embryo. The results are shown 
in Fig.[2t: Hb boundary precision is lost, and many pro- 
files shift towards the anterior Jxjrb = 0.43 (±0.05) EL), 
in agreement with experiment j3(. 

To isolate which emhrvo-to-emhrvo variations are most 



0.12 
0.1 
0.08 
£ 0.06 
0.04 
0.02 




g-g/ b and/ R 


uncorrelated 




b-h J B and J R 


correlated 




- l*Hb 






■ ^^hb 






a 












K 













0.1 0.2 0.3 0.4 0.5 
relative magnitude of fluctuations 

FIG. 3: Magnitude of variations in position of Hb bound- 
ary as a function of relative magnitude (standard deviation 
a/mean value) of embryo-to-embryo variations in model pa- 
rameters. Both uncorrelated and correlated variations in the 
injection rates Jb, Jr are shown (see text). When varying the 
fluctuations of one parameter, other parameters were fixed to 
their wild type values, indicated by closed symbols. 

important, we studied the magnitude of the fluctuations 



in the Hb boundary position as a function of the embryo- 
to-embryo variability of various parameters. Correlated 
variations that affect Bed and the co-repressor similarly 
have little effect, due to the symmetry of the activator- 
co-repressor geometry (see Fig. EJ. If, for instance, both 
Bed and co-repressor are degraded by the same pathway, 
then the embryo-to-embryo variations in their degrada- 
tion rates are correlated and thus filtered out (see Fig. [21 . 
Furthermore, it is conceivable that not only degrada- 
tion, but also synthesis of Bed and co-repressor exhibit 
correlated embryo-to-embryo fluctuations, because they 
share, for instance, the same components of the tran- 
scription/translation apparatus. These correlated vari- 
ations of the Bcd/co-repressor injection rates are again 
filtered (see Fig. [3J . However, Bcd/co-repressor synthesis 
may also exhibit some uncorrelated variation. As we dis- 
cussed above, our symmetric activator-co-repressor sys- 
tem buffers fluctuations far more than would be possible 
in a single morphogen system. Nevertheless, these uncor- 
related fluctuations are not completely filtered out and, if 
large enough, can adversely affect the accurate position- 
ing of the hb expression domain (see Fig.[3|. Our model 
therefore predicts that Bed synthesis must be tightly con- 
trolled, at least relative to that of the co-repressor. 
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